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Preface 


This  report  describes  the  computer  program  to  produce  the  terrain  data  for 
the  Battlescale  Forecast  Model  (BFM)  using  the  U.  S.  Geological  Survey's 
(USGS)  30  arc-second  elevation  data  (GTOPO30).  The  program  written  in 
the  C  language  is  also  attached  in  appendix  A. 
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Executive  Summary 

A  new  computer  program  named  read_terrain.c  is  developed  to  generate 
the  elevation  data  for  the  Battlescale  Forecast  Model  (BFM)  from  the  U.  S. 
Geological  Survey  (USGS)  30  arc-second  elevation  data  (GTOPO30  data  set). 
Since  the  data  set  covers  the  entire  globe,  terrain  data  for  the  BFM  can  be 
obtained  for  any  part  of  the  world,  except  for  the  areas  in  high  latitudes. 

This  report  describes  the  GTOPO30  data  set,  which  contains  the  data  every 
30  arc-second  over  the  entire  globe.  The  data  was  originally  contained  in 
five  compact  disk  read  only  memory  (CDROM)s  to  increase  the  portability 
and  ease  of  use.  The  data  were  compressed  and  stored  on  a  single  CDROM. 

The  method  of  data  extraction  from  the  data  set  is  mentioned,  which  is 
capable  of  producing  the  terrain  data  for  the  BFM  from  a  multiple  number 
(up  to  four)  tiles.  Several  function  subprograms  are  used  in  this  "C" 
language  program,  and  the  role  of  each  subprogram  is  described.  Examples 
of  terrain  data  produced  by  this  program  are  shown.  The  program  codes  in 
"C"  language  are  included  in  appendix  A. 
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1.  Introduction 


Battlescale  Forecast  Model  (BFM),  developed  at  the  U.S.  Army  Research 
Laboratory  (ARL),  has  been  extensively  used  to  produce  short-range 
forecasts  of  atmospheric  conditions  as  a  component  in  both  the  Integrated 
Meteorological  System  (IMETS)  and  the  Computer  Assisted  Artillery 
Meteorology  (CAAM)  system. 

1.1  Background 

The  elevation  data  is  an  essential  input  data  for  a  mesoscale  weather 
forecast  model  such  as  the  BFM.  Previously,  the  BFM  has  obtained  elevation 
data  from  the  Digital  Terrain  Elevation  Data  (DTED)  Level  1  data  set, 
produced  by  the  National  Imagery  and  Mapping  Agency  (formerly  the 
Defense  Mapping  Agency).  DTED  Level  1  has  a  horizontal  grid  spacing  of 
3-arc  seconds  (approximately  90  m).  The  data  on  the  compact  disk  read 
only  memory  (CDROM)  usually  covers  an  area  of  12°  of  latitude  by  18°  of 
longitude.  Therefore,  a  substantial  number  of  CDROMs  are  needed  to  cover 
the  entire  globe.  The  BFM  is  typically  used  over  model  areas  ranging  from 
100  X  100  km  to  500  x  500  km  with  a  model  grid  spacing  of  5  to  10  km.  For 
the  current  BFM  on  IMETS,  elevation  data  for  an  area  of  1,600  x  1,600  km 
is  needed  for  a  model  domain  of  500  x  500  km.  Therefore,  depending  on  the 
area  of  interest,  a  number  of  CDROMs  are  required  to  create  elevation  data 
for  the  model.  A  serious  shortcoming  of  the  DTED  Level  1  data  set  is  that 
there  are  regions  where  no  data  are  provided,  most  notably  various  regions 
in  the  Southern  Hemisphere. 

The  U.S.  Geological  Survey  (USGS)  GTOPO30  is  a  global  digital  elevation 
data  set  spanning  the  entire  globe.  The  elevation  data  in  GTOPO30  are 
regularly  spaced  at  30-arc  seconds  (approximately  1  km).  Thus,  it  is  an 
appropriate  data  source  to  create  elevation  data  for  BFM  application  in  any 
part  of  the  world. 

1.2  Purpose 

The  purpose  of  this  technical  note  is  to  describe  the  method  and  a  computer 
program  developed  to  create  the  elevation  data  for  BFM  from  the  GTOPO30 
data  set: 

•  section  H,  GTOPO30  data  are  described, 

•  section  III,  the  method  of  data  reduction  is  described, 

•  section  FV,  the  role  of  the  function  subprogram  in  the  program  is 
described, 

•  section  V  examples  of  terrain  data  produced  by  the  program  are  shown, 

•  section  VI  briefly  mentions  limitations  of  the  program, 

•  section  VII  is  a  summary,  and 

•  appendix  A,  the  computer  program  in  C  language  is  provided. 
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2.  GTOPO30  Data  Set 


GTOPO30  is  a  global  data  set  covering  the  full  latitude  range  (90°  south  to 
90°  north),  and  the  full  range  of  longitude  (180°  west  to  180°  east).  The 
horizontal  grid  spacing  is  30-arc  seconds.  The  data  is  given  in  meters  above 
mean  sea  level  (MSL).  Ocean  areas  are  assigned  a  value  of  -9999.  The  data 
is  provided  as  16-bit  signed  integer  data  in  a  simple  binary  format. 

GTOPO30  is  divided  into  33  small  pieces  or  tiles.  The  area  from  60°  south 
latitude  to  90°  north  latitude  and  from  180°  west  longitude  to  180°  east 
longitude  is  covered  by  27  tiles,  with  each  tile  covering  50°  of  latitude  and 
40°  of  longitude. 

Antarctica  (90°  south  latitude  to  60°  south  latitude  and  180°  west  longitude 
to  180°  east  longitude)  is  covered  by  6  tiles.  Each  tile  covers  30°  of  latitude 
and  60°  of  longitude.  Figure  1  shows  the  distribution  of  tiles  covering  the 
entire  globe. 

Each  of  the  27  tiles  that  individually  cover  50°  of  latitude  and  40°  of 
longitude  have  6,000  rows  and  4,800  columns.  Each  of  the  six  Antarctica 
tiles  that  individually  cover  30°  latitude  and  60°  of  longitude  have  3,600 
rows  and  7,200  columns. 

Details  of  GTOPO30  are  described  in  the  documentation  found  in  the 
universal  resource  locator  (URL)  address  on  the  web: 

http:/ / edcAvww.cr.usgs.gOv/landdaac/gtopo30/ 

To  increase  the  portability  and  ease  of  use  of  this  software,  the  data  for  33 
tiles  have  been  extracted  from  5  CDROMs  obtained  from  the  USGS  and 
compressed  using  the  gzip  command  enabling  the  data  for  33  tiles  to  be 
stored  on  a  single  CDROM.  Through  the  use  of  the  gzip  utility,  the  total 
size  of  the  global  data  is  reduced  from  almost  2.72  GB  to  approximately 
290  MB.  Thus,  elevation  data  for  any  part  of  the  world  can  be  obtained  from 
a  single  CDROM.  Table  1  gives  the  names  of  the  tiles  contained  in  the 
CDROM  and  the  USGS  names  of  the  tiles. 
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Figure  1.  Map 
showing  GTOPO30 
tiles  (obtained  from 
the  GTOPO30  home 
web  page). 


Table  1.  Tile  and  USGS  names 


Number 

Tile 

name 

USGS 

names 

1 

e020n9~l 

E020N90.DEM 

2 

e060n9~l 

E060N90.DEM 

3 

el00n9~l 

E100N90.DEM 

4 

el40n9~l 

E140N90.DEM 

5 

wl80n9~l 

W180N90.DEM 

6 

wl40n9~l 

W140N90.DEM 

7 

wl00n9~l 

W100N90.DEM 

8 

w060n9-'l 

W060N90.DEM 

9 

w020n9~l 

W020N90.DEM 

10 

e020n4~l 

E020N40.DEM 

11 

e060n4~l 

E060N40.DEM 

12 

el00n4~l 

E100N40.DEM 

13 

el40n4~l 

E140N40.DEM 
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Table  1.  Tile  and  USGS  names  (continued) 


Number 

Tile 

name 

USGS 

names 

14 

wl80n4~l 

W180N40.DEM 

15 

wl40n4~l 

W140N40.DEM 

16 

wl00n4~l 

W100N40.DEM 

17 

w060n4~l 

W60N40.DEM 

18 

w020n4~l 

W020N40.DEM 

19 

e020sl~l 

E020N40.DEM 

20 

e060sl~l 

E060S10.DEM 

21 

el00sl~l 

ElOOSlO.DEM 

22 

el40sl~l 

E40S10.DEM 

23 

wl80sl~l 

W180S10.DEM 

24 

wl40sl~l 

W140S10.DEM 

26 

w060sl~l 

W060S10.DEM 

27 

w020sl~l 

W020S10.DEM 

28 

e060s6~l 

E060S60.DEM 

29 

el20s6~l 

E120S60.DEM 

30 

wl80s6~l 

W180S60.DEM 

31 

wl20s6~l 

W120S60.DEM 

32 

w060s6~l 

W060S60.DEM 

33 

w000s6~l 

W000S60.DEM 

In  table  1,  each  name  represents  the  latitude  and  longitude  of  the  northwest 
comer  of  tiles.  The  names  of  the  USGS  files  changed  as  they  were  gzipped 
and  stored  on  CDROM  due  to  file  name  length  limitations  in  the  disk 
operating  system  (DOS).  For  instance,  tile  #  7,  wl00n9~l  covers  longitudes 
from  100°  west  to  60°  west,  and  latitudes  from  90°  north  to  40°  north. 


3.  Method  of  Elevation  Data  Creation  for  BFM 


The  procedures  to  create  elevation  data  sets  from  the  CDROM  are  as 
follows: 

In  the  BFM,  Cartesian  coordinates  are  used  for  horizontal  coordinates. 
Therefore,  the  latitude  and  longitude  of  model  grid  points  (I,  J)  are 
calculated  as 


e(j)  =  eo+ 


j- 


(1) 


^(l,  j)  =  (j)0  + 


A- 

'imox  +  lV 
2  ). 

1 

)  -  cos 

((180)®^]] 

(2) 


where 

R. 


4>o 

eo 


(linax5  ^max) 

(1,J) 


the  radius  of  the  earth 
the  unit  grid  distance 

longitude  of  the  center  of  a  model  domain 
latitude  of  the  center  of  model  domain 

the  mamimum  numbers  of  grid  points  for  the  (x,  y)  directions 
grid  point  counters  in  the  x  and  y  directions 


Therefore,  minimum  and  maximum  longitudes  and  latitudes  of  the  model 
domain  are  given  as 

For  00^0, 


minimum  longitude  =  (})„i„  =  4)(1,  J„„) 
maximum  longitude  =  =4>(linax,  Jmax) 

and  for  6o<0, 


minimum  longitude  =  (j)mi„=  1) 

maximum  longitude  =  (j)„„=  1) 

minimum  latitude  =  =  0(1,  1)  =  0(lmax>  1) 
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maximum  latitude  =  =  0(1,  J„ax)  =  0(lmax>  ^max) 

Since  the  data  are  given  at  every  30  arc  second  (or  1/120°),  the  numbers  of 
a  data  point  in  longitudinal  and  latitudinal  directions,  i,  and  jt  for  an  area 
bounded  by  4)„^,  0„^,  and  0„3,  are  calculated  as 


it  =  (4>max-4)min)x  120+1 


j.=  (emax-0mJx  120+1 


(3) 

(4) 


The  longitudes  and  latitudes  of  the  elevation  data  covering  the  model 
domain  are  given  by 


4)(i)  =  <l)n,in+  (i  -1)  X  (1/120),  i  =1,  2, . ,i,  (5) 

0(j)  =  e.i„  +  (i-l)x(l/12O),  j=l,2, . ,j,  (6) 

All  tiles  are  examined  to  see  if  they  contain  any  elevation  data  for  (j)(  i ), 

i=l,2 . ,  i, ,  and  0(  j ),  j=l, . ,  j,  and,  any  tile(s)  containing  the  data  are 

copied  to  the  directory  in  which  the  program  is  being  executed  and 
uncompressed.  Depending  on  the  location  and  size  of  the  model  domain, 
elevation  data  may  have  to  be  obtained  from  a  multiple  number  of  tiles. 

The  present  program  is  designed  so  that  elevation  data  can  be  obtained  from 
1,  2,  or  4  tiles.  Therefore,  if  the  model  domain  is  located  at  high  latitudes 
where  longitudmal  lines  are  densely  packed,  the  present  program  may 
become  incapable  of  extracting  elevation  data,  depending  on  the  area  size. 
The  limitations  of  the  program  will  be  discussed  in  sectionVI. 

Within  each  tile,  the  data  is  stored  in  row  major  order  (all  the  data  for  row 
1,  followed  by  the  data  for  row  2,  etc.).  The  northwest  comer  is  designated 
as  (0, 0)  in  the  program.  The  portion  of  data  in  the  tiles  which  encompass  the 
model  domain  are  extracted  and  used  to  interpolate  elevation  data  for  the 
model  domains. 

If  a  grid-point  (x',  y')  of  the  model  domain  is  siurrormded  by  four  grid 
points  of  the  elevation  data,  an  interpolated  elevation  value  Z'  at  (x',  y')  is 
then  calculated  using  bilinear  interpolation  as 

t,  =  Z(x,  y)  +  (X'  -  x)-[  Z(x  +  1 ,  y)  -  Z(x,  y)]  (7) 

t^  =  Z(x,  y  +  1)  +  (x'  -  x)-[  Z(x  +  1,  y  +  1)  -  Z(x,  y  +  1)]  (8) 

Z'(x',y')  =  t,  +  (y'-y)-[l2-t,]  (9) 

Here,  (x,  y)  is  the  southwest  grid  of  four  grid  points  surrounding  the  model 
grid  point  (x',  y '),  and  Z  (x,  y)  is  the  elevation  data  at  (x,  y).  This  operation 
is  repeated  for  all  grid  points  of  the  model  domain. 


4.  Program  Description 

Several  functions  subprograms  are  used  in  this  program.  In  the  following, 
the  role  of  each  function  is  described: 

•  main  The  latitude  and  longitude  of  the  center  of  the  model  domain,  the 
number  of  grid  points  in  the  x  and  y  directions,  and  the  grid  spaces  in 
kilometers  are  typed  in  this  program.  Latitude  and  longitude  values  of  the 
model  domain  grid  points  are  calculated,  and  maximum  and  minimum 
latitudes  and  longitudes  are  passed  to  the  function  called  dataSOs. 

•  dataSOs  This  subprogram  examines  which  and  how  many  tile(s)  encompass 
the  model  domain.  The  tiles  containing  necessary  data  are  copied,  and 
uncompressed  by  the  gzip  -d  command  to  the  directory  where  the  program 
is  executed.  Elevation  data  on  the  tile(s)  are  read  from  the  northwest  comer; 
(for  example,  all  the  data  for  row  1,  followed  by  all  the  data  for  row  2).  After 
one  row  of  the  data  is  read,  the  data  is  then  passed  onto  the  fvmction 
get_datal  or  get_data2  in  order  to  select  the  portion  of  data  that 
encompasses  the  model  domain. 

•  get_datal  and  get_data2  These  functions  examine  whether  the  elevation 
data  encompasses  the  model  domain.  An  array  of  elevation  data  is  then 
made  in  these  functions  and  passed  onto  the  function  zinterp.  The  function 
get_datal  is  used  for  the  tiles  1  to  27  and  get_data2,  for  the  tiles  28  to  33. 

•  zinterp  Is  a  subprogram  that  performs  bilinear  interpolation  of  elevation 
data  for  the  grid  points  of  the  model  domain  and  creates  "terrain_data"  for 
the  model  dornain. 

•  cp_from_cd  The  data  on  a  tile  are  copied  from  a  CDROM  onto  the  directory 
where  the  program  resides. 

•  unzip  To  unzip  the  command  gzip  -d  is  executed  to  uncompress  data  for 
a  given  tile. 

•  rm_file  After  uncompressed  data  have  been  read,  data  file(s)  are  removed 
from  the  directory  of  program  execution. 
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5.  Examples  of  Terrain  Data 

The  computer  program  is  capable  of  producing  terrain  data  for  the  BFM 
from  some  multiple  numbers  of  tiles.  In  this  section,  three  different  terrain 
data  files  produced  from  different  numbers  of  tiles  are  presented; 

•  Western  Alaska  and  Bering  Strait  Area  Figure  2  shows  the  terrain  data  for 
the  area  of  1,000  x  1,000  km  (101  x  101  grids)  centered  at  65.0°  N  and  165.0° 
W,  and  grid  spacing  of  10  km.  These  data  are  obtained  from  one  tile, 
wl80n9~l.  The  DTED  level  1  does  not  have  complete  coverage  of  this  area. 

•  Brazil  Area  The  DTED  Level  1  data  set  does  not  have  any  data  coverage 
over  Brazil.  By  using  the  present  program  and  data  set,  the  terrain  data  over 
Brazil  can  be  generated.  Figure  3  represents  the  central  part  of  Brazil,  the 
area  centered  at  8.0°  S  and  39.5°  W  covering  1,000  x  1,000  km  (51  x  51  grids) 
with  grid  spacing  of  20  km. 

•  Midwestern  US  Figure  4  shows  the  terrain  data  obtained  from  four  tiles, 
wl40n9~l,  wl00n9~l,  wl40n4~l,  and  wl00n4~l.  The  center  of  this  domain 
is  located  at  the  four  comer  point,  40°  N  and  100°  W.  The  data  covers 
1,000  km  X  1,000  km  (101  x  101  grids)  with  10  km  grid  spacing. 
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Figure  2.  Terrain  data  showing 
northwestern  Alaska  and  the 
eastern  edge  of  Siberia,  centered 
at  65°  N  and  165°  W.  Data  were 
obtained  from  one  tile. 


CONTXtUR  FROM  100  TO  1700  BV  100 


Figure  3.  Terrain  data  for 
eastern  Brazil  entered  at  8°  S  and 
39.5°  W.  Data  were  obtained 
from  two  tiles. 
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Figure  4.  Terrain  data  for  the 
midwestem  United  States  centered 
at  40.0°  N  and  100°.  Data  were 
obtained  from  four  tiles. 


6.  Limitation  of  the  Program 

Because  the  program  is  capable  of  extracting  elevation  data  only  from  two 
consecutive  tiles  in  both  longitudinal  and  latitudinal  directions,  data  in  high 
latitudes  cannot  be  obtained. 

Table  2  shows  the  limits  of  latitudes  over  which  elevation  data  cannot 
be  obtained. 


Table  2.  Latitude  limits 


In  the  attached  program  code  (appendix  A)  the  greatest  number  of  grids  for 
the  model  domain  is  currently  set  to  161  x  161.  Therefore,  if  this  grid 
number  is  exceeded,  the  program  must  be  modified.  The  program  is 
developed  for  the  UNIX  operating  system.  One  additional  limitation  is  that 
the  terrain  data  with  grid  spacing  smaller  than  30  arc-second  (about 
1  km)  cannot  be  produced. 
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7.  Summary 


The  present  program  is  developed  to  produce  terrain  data  for  the  BFM  from 
the  uses  GTOPO30  (global  digital  elevation  data  set).  GTOPO30  data  are 
extracted  from  five  CDROMs,  compressed  by  the  UNIX  utility  "grip"  and 
archived  onto  a  single  CDROM.  Thus,  terrain  data  for  most  of  the 
BFM  application  can  be  obtained  from  the  CDROM  by  using  the 
present  program. 
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Acronyms 


ARL 

U.S.  Army  Research  Laboratory 

BFM 

Battlescale  Forecast  Model 

CAAM 

Computer  Assisted  Artillery  Meterology 

CDROM 

compact  disk  read  only  memory 

DOS 

disk  operating  system 

DIED 

Digital  Terrain  Elevation  Data 

IMETS 

Integrated  Meterological  System 

MSL 

mean  sea  level 

URL 

universal  resource  locator 

ARL 

U.S.  Army  Research  Laboratory 

uses 

U.S.  Geological  Survey 

Appendix  A 

Program:  read_terrain.c 
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/*  read  terrain.c  */ 

/*  Program  to  extract  and  interpolate  an  elevation  */ 

/*  data  file  from  up  to  4  USGS  GTOPO30  data  tiles,  depending  */ 
/*  on  how  many  boundaries  are  overlapped.  In  *! 

I*  these  files,  -999  represents  water  so  it  is  */ 

/*  replaced  with  0  as  it  is  read.  */ 

/*  Steve  Kirby  (30  July  1999)  */ 

#include  <stdio.h> 

#include  <math.h> 


#define  ip  161 
#define  jp  161 
#define  i_x  4000 
#define  j_y  4000 
#define  nl  4800 
#define  n2  6000 
#define  n3  7200 
#define  n4  3600 

#define  conv_factor  3 . 1 4 1 6/1 80 . 

#define  re  6371.22 

FILE*fp,*fpl,*fp2,*fp3,*fp4,*fpout,*fpcheck; 
float  longi,lati,gridsp; 
float  origlon; 

int  val,min_val,max_val,numx,numy; 
int  i,j,imax=4800,jmax=6000,ipl,jpl; 
float  lat_b[jp],lon_b[ip][jp],gdis,zs[ip][jp],lato,lono; 
float  center=l. 0/120.0; 

float  z[i_x][j_y],lat[j_y],lon[i_x],hi_2[i_x][j_y], 
lat_o_2[j_y]; 

float  Ion_o_2[i_x],lon_o2[2*n3]; 
signed  short  hi[2*nl]; 
signed  short  hi2[2’^n3]; 
int  i_lon,j_lat; 

float  xminlon,xmaxlon,xminlat,xmaxlat; 
float  min_lat,min_lon,max_lat,max_lon; 
float  lon_ol[2*nl]; 

/*  protot5Yes  *  / 

void  data30s(float,float,float,float,float); 
void  get_datal  (float,float,float,float,int,int,float. 


float,mt*); 

void  get_data2(float,float,float,float,int,int,float 
float,mt*); 

void  zinterp(int,int,float  [],float  []); 
void  cp_from_cd(char[]); 
void  unzip(char  []); 
void  nn_file(char  []); 
int  main(){ 

/*  user  input  for  center  of  domain  */ 
printfC'enter  the  center  lat/lon  of  the  domain 
(-inSHandWH):"); 
scanf("%f  %f",&lato,&lono); 

I*  need  to  save  input  Ion  for  output  */ 

/*  in  case  it's  negative  * ! 
origlon=lono; 

/*  user  input  for  size  of  grid  * ! 
printfC'enter  the  number  of  grid  points  in  the 
X-  and  y-direction: "); 
scanf("%d  %d",&ipl,&jpl); 

/*  user  input  for  grid  spacing  in  kms  * ! 
printfC’enter  the  grid  spacing  (kms): "); 
scanfC’%f",&gdis); 

printfC'\n\ncenter  lat  I  lon:%5.2f  %6.2f\n 

size  of  grid:%dx%d\n  grid  spacing:%6.2f  kms\n",lato,lono4pl,jpl,gdis) 

I*  formulate  lat  1  Ions  of  BFM  grid  points  */ 
for(j=0;j<jpl;j++){ 

lat_b[j]=lato  +  ((float)(2*(j+l)-jpl-l)/2.)  *  (gdis/(re*conv_factor)); 
for(i=0;i<ipl;i++){ 

lon_b[i]  [j]=lono+(float)(2*(i+l)-ipl-l)  /2.*(gdis  /(re*convJactor* 
cos(lat_b[j]*conv_factor))); 

} 

} 

min_lat=lat_b[0]; 

max_lat=lat_b[jpl"l]; 

if(lato<0.0) 

{ 

max_lon=lon_b[ipl-l]  [0]; 
min_lon=lon_b  [0]  [0]; 

} 

else 

{ 

min_lon=lon_b[0][jpl-l]; 

max_lon=lon_b[ipl-l][jpl-l]; 

} 


data30s(max_lat,niin_lat,max_lon,min_lon,center); 
return  0; 

} 

void  data30s(float  xmaxlat, float  xminlat,  float 
xmaxlon, float  xminlon, float  center) 

{ 

float  dlat0[33],d]on0[33],dlatl[33],dlonl[33]; 

signed  short  hi_max=-1000,hi_min=1000; 

int  kd[33]; 

int  lonjndex; 

int  nrsdt,jt,nk,k,i,j; 

float  lon_f,lat_f; 

char  lon_f_str[4],lat_f_str[3]; 

float  lat[j_y],lon[i_x]; 

float  xlonO,xlatO,xlonl,xlatl; 

float  XX, yy; 

int  imaxl,jmaxl,imax2,jmax2,imax,jmax; 
float  lat_ol[2*n3],lat_o2[2*n3]; 
intindex_region,index_r,index,i_switch,j_switch; 
float  dlat_max,dlat_min,dlon_max,dlon_min; 
char  *elev_data[4]; 
intjj; 

float  dlon00[4]; 
char  *dem_file[]={ 

"e020n9~l","e060n9~l","el00n9~l", 

"el40n9~l”, "wl80n9~l", "wl40n9~l", 

"Wl00n9~l","w060n9~l","w020n9~l", 

"e020n4~l",”e060n4~l","el00n4~l", 

"el40n4~l","wl80n4~l","wl40n4~l”, 

"Wl00n4~l","w060n4~l","w020n4~l", 

"e020sl~l","e060sl~r',"el00sl~l", 

"el40sl~l","wl80sl~l","wl40sl~l", 

"wl00sl~l","w060sl~l","w020sl~l", 

"e060s6~l","el20s6~l","wl80s6~l", 

"wl20s6~l","w060s6~l","w000s6~l"}; 

char  show_file[60]; 

printf("\nbounding  region--  LAT:%5.2f  %5.2f  LON:%5.2f 

%5.2f\n",xininlat,xmaxlat, 

xininlon,xmaxlon); 

nrs=(int)(l. /center  +  0.5); 

it=(xinaxlon-xminlon)*nrs; 

jt=(xmaxlat-xminlat)*nrs; 


nk=-l; 

for(k=0;k<33;k++){ 
lon_f_str[0]=dem_file[k]  [1]; 
lon_f_str[l]=dem_file[k][2]; 
lon_f_str[2]=dem_file[k]  [3]; 
lon_f_str[3]='\0'; 
lon_f=(float)atoi(lon_f_str); 
lat_f_str[0]=dem_file[k]  [5]; 
lat_f_str[l]='0’; 
latj_str[2]='\0'; 
lat_f=(float)atoi(lat_f_str); 
if(dem_file[k][0]  ==  'w') 
xlonO=-lon_f; 
else 

xlonO=lon_f; 

if(dem_file[k][4]  ==  ’s') 
xlatl=-lat_f; 

else 
xlatl=lat_f; 
if(k  <=  26){ 

xlonl=xlon0+(float)(nl-l)/120.; 
if(xlonl  >=  360.0) 
xlonl-=360.; 

xlat0=xlatl-(float)(n2-l)/120.; 

} 

else! 

xlonl=xlon0+(float)(n3-l)/120.; 
if(xlonl  >=  360.0) 
xlonl-=360.; 

xlat0=xlatl-(float)(n4-l)  / 120.; 

) 

if((xinirLlon  >  xlonl)  &&  (xminlat  >  xlatl)) 
goto  label20; 

if((xminlon  ==  -180.)  &&  (xmaxlon  ==  180.)){ 

I*  the  pole  is  inside  the  domain  * ! 

if((xmaxlat  ==  90.)  &&  (k  <=  8)){ 

nk+=l; 

kd[nk]=k; 

dlat0[nk]=xlat0; 

dlon0[nk]=xlon0; 

dlatl[nk]=xlatl; 

dlonl  [nk] =xlonl; 

} 

else  if((xminlat  ==  -90.)  && 

(k>26)){ 

/*  South  Pole  V 


nk+=l; 

kd[nk]=k; 

dlatO[nk]=xlatO; 

dlonO[nk]=xlonO; 

dlatl[nk]=xlatl; 

dlonl[nk]=xlonl; 

) 

) 

else! 

I*  the  pole  is  outside  the  domain  */ 

for(i=0;i<it;i++){ 

xx=xminlon+i*center; 

for(j=0;j<jt,j++){ 

yy=xminlat+j*center; 

if((xx>=xlonO)  && 

(xx<xlonl)  && 

(yy>xlatO)  && 

(yy<xlatl)  && 

(y>=xminlat)  && 

(yy  <  xmaxlat)  && 
rijc+=l; 

dlatO[nk]=xlatO; 

dlonO[nk]=xlonO; 

dlatl[nk]=xlatl; 

dlonl[nk]=xlonl; 

kd[nk]=k; 

goto  label20; 

) 

) 

) 

) 

label20: 

continue; 

) 

printf("\nneed  %d  30-sec  tile(s):\n",nk+l); 

for(k=0;k<=nk;k++){ 

strcpy(show_file,dem_file[kd[k]]); 

show_file[6]='0'; 

show_f ile  [7] = '  \  O'; 

printf("%d»  %s  LAT:  %5.2f  %5.2f  LON: 
%6.2f  %6.2f\n",k+l,show_file,dlat0[k], 
dlatl  [k],dlonO[k],dlonl  [k]); 

} 
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lon_index=l; 

if((nk  ==  1  &&  dlon0[0]  ==  20.0  &&  dlon0[l] 

==  -20.0)  I  I  (nk  ==  1  &&  dlon0[0]  ==  60.0  &&  dlon0[l]  ==  0.0)){ 
lon_index=2; 

elev_data[0]=dem_file[kd[l]]; 

elev_data[l]=dem_file[kd[0]]; 

dlon00[0]=dlon0[l]; 

dlon00[l]=dlon0[0]; 

if(nk  ==2  &&  dlon0[l]  ==  60.0  &&  dlon0[2] 

==0.0)( 

lon_index=2; 

elev_data[0]=dem_file[kd[0]]; 

elev_data[l]=dem_file[kd[2]]; 

elev_data[2]=dem_file[kd[l]]; 

dlon00[l]=dlon0[2]; 

dlon00[2]=dlon0[l]; 

} 

else  if(nk  ==  3  &&  dlon0[0]  ==  20.0  && 
dlon0[l]  ==  -20.0  &&  dlon0[2]==  20.0  && 
dlon0[3]==-20.0){ 
lon_index=2; 

elev_data[0]=dem_file[kd[l]]; 

elev_data[l]=dem_file[kd[0]]; 

elev_data[2]=dem_file[kd[3]]; 

elev_data[3]=dem_file[kd[2]]; 

dlon00[0]=dlon0[l]; 

dlon00[l]=dlon0[0]; 

) 

else  if(nk  ==  3  &&  dlon0[2]  ==  60.0 
&&  dlon0[3]  ==  0.0)( 
lon_index=2; 

elev_data[2]=dem_file[kd[3]]; 

elev_data[3]=dem_file[kd[2]]; 

dlon00[2]=dlon0[3]; 

dlon00[3]=dlon0[2]; 

) 

else{ 

for(k=0;k<=nk;k++) 

elev_data[k]=dem_file[kd[k]]; 

) 

} 

else{ 
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for(k=0;k<=nk;k++) 

elev_data[k]=dem_file[kd[k]]; 

) 

for(k=0;k<=nk;k++){ 
if(dlatO[k]  <  -60.){ 
index_region=2; 
goto  label22; 

} 

else( 

index_region=l; 

jmax=n2; 

imax=nl; 

} 

) 

label22: 

if(index_region  ==  1){ 
if(nk  ==  0){ 

I*  lat  I  Ion  of  elevation  data  */ 

labelSOOO; 

for(j=0;j<jmax;j++) 

lat_ol[j]=dlatO[nk]+ 

center’^(float)j; 

for(i=0;i<imax;i++) 

lon_ol  [i] =dlonO  [nk] + 

center’^(float)i; 

/*  open  the  elevation  data  *  j 

cp_from_cd(elev_data[0]); 

unzip(elev_data[0]); 

if  ( !  (fp=f  open(elev_data  [0],"r")))  { 

printfC'USGS  DEM  file 

%s  unavailable ... 

exiting. .. \n",dem_file  [kd  [nk]  ] 

); 

exit(O); 

} 

j_lat=-l; 

j=jmax-l; 

while(j>=l)| 

for(i=0;i<innax;i++){ 

fread(&hi[i], 

sizeof(signed  short),l,fp); 
if(hi[i]  ==  -9999) 
hi[i]=0; 

if(hi[i]>hi_max) 

hi_max=hi[i]; 

if(hi[i]<hi_min) 

hi_min=hi[i]; 


index_r=l; 

get_datal  (xminlat,xmaxlat, 

xminlon,xmaxlon,0,iniaxdat_ol  [j],lat_ol  [j-l],&index_r); 
if(index_r  ==  2) 
goto  label210; 

j--; 

) 

label210: 

fclose(fp); 

nn_file(elev_data[0]); 

) 

else  if(nk  ==  1){ 

I*  tile  1  is  above,  and 
tile  2  is  below  *  j 
if(dlon0[0]  ==  dlonO[l]  && 
dlonl[0]  ==  dlonl[l]){ 
printfC'case  tile  1  above, 
tile  2  belowXn"); 
f  or(j=0;j  <2*jmax;j ++) 
lat_ol[j]=dlatO[l]+ 
center*(float)j; 

,  for(i=0;i<imax;i++) 
lon_ol  [i]=dlonO[l]+ 
center*(float)i; 

/*  open  elevation  data  */ 

cp_from_cd(elev_data[0]); 

unzip(elev_data[0]); 

cp_from_cd(elev_data[l]); 

unzip(elev_data[l]); 

if(!(fpl= 

fopen(elev_data[0],"r"))){ 

printf("USGS  DEM  file  %s  unavailable  ...  exiting  ...\n",elev_data[0]); 
exit(O); 

} 

if(!(fp2= 

fopen(elev_data[l],"r"))){ 

printffUSGS  DEM  file  %s  unavailable  ...  exiting  ...\n",elev_data[l]); 
exit(O); 

} 

j_lat=-l; 

j=2*jmax-l; 

while(j>=l){ 

if(j>=jmax){ 

for(i=0;i<imax;i++){ 

fread(&hi[i]. 


sizeof(signed  short),l,fpl); 
if(hi[i]  ==  -9999) 
hi[i]=0; 


) 

else{ 

for(i=0;i<imax;i++){ 

fread(&hi[i], 

sizeof(signed  short),l/fp2); 
if(hi[i]  == 

-9999) 

hi[i]=0; 

} 

) 

index_r=l; 

get_datal(xminlat,xniaxlat, 

xminlon,xmaxlon,0,imax, 

lat_ol[j], 

lat_ol  [j-1  ],&index_r); 
if(index_r  ==  2) 
goto  labelSlO; 

j-; 

} 

labelSlO: 

fclose(fpl); 

rm_file(elev_data[0]); 

fclose(fp2); 

rm_file(elev_data[l]); 

} 

I*  tile  1  is  to  the  left; 
tile  2  is  to  the  right  * ! 
else  if(dlat0[0]  ==  dlatO[l] 
&&  dlatl[0]  ==  dlatl[l]){ 
printfC'case  tile  1  to 
left  tile  2  to  rightXn"); 
labelbOOO: 
for(j=0;j<jmax;j++) 
lat_ol[j]=dlat0[0] 

+  center*(float)j; 

if(lon_index  ==!){ 
for(i=0;i<2*imax;i++) 
lon_ol[i]= 
dlon0[0]  + 
center*(float)i; 

} 

else  if(lon_index  ==  2){ 
for(i=0;i<2*imax;i++) 


lon_ol[i]= 

dlon00[0]  +  center*(float)i; 


I*  open  the  elevation  data  */ 
cp_from_cd(elev_data[0]); 
unzip(elev_data[0]); 
cp_from_cd(elev_data  [  1  ]); 
unzip(elev_data[l]); 

if(!(fpl= 

fopen(elev_data[0],"r")))| 

printf("USGS  DEM  file  %s  unavailable  ...  exiting  ...\n",elev_data[0]); 
exit(O); 


if(!(fp2= 

fopen(elev_data[l],"r"))){ 

print  fC’USGS  DEM  file  %s  unavailable  ...  exiting  ...\n",elev_data[l]); 
exit(O); 

) 

j_lat=-l; 

j=jmax-l; 

while(j>=l){ 

for(i=0;i<imax;i++){ 

fread(&hi[i], 

sizeoffsigned  short),l/fpl); 
if(hi[i]  ==  -9999) 
hi[i]=0; 

) 

for(i=imax;i<2*imax;i++){ 

fread(&hi[i], 

sizeoffsigned  short),l,fp2); 
if(hi[i]  ==  -9999) 
hi[i]=0; 

) 

index_r=l; 

get_datal(xminlat,xmaxlat, 

xminlon,xmaxlon,0,2*imax, 

lat_ol[j], 

lat_ol  [j-1  ],&index_r); 
if(index_r  ==  2) 
goto  label410; 


Iabel410: 

fclose(fpl); 

rm_file(elev_data[0]); 

fclose(fp2); 

rm_file(elev_data[l]); 

} 

} 

else  if(nk  ==  3){ 

/*  4-tiles  * ! 

for(j=0;j<2*jmax;j++) 

lat_Dl[j]=dlatO[2]+ 

center*(float)j; 

if(lon_index  ==  1){ 

for(i=0;i<2*imax;i++) 

lon_ol[i]=dlon0[0] 

+  center*(float)i; 

} 

else  if(lon_index  ==  2){ 
for(i=0;i<2*imax;i++) 
lon_ol[i]= 
dlon00[0]  + 
center  *  (float)i; 

} 

I*  open  elevation  data  * ! 
cp_from_cd(elev_data[0]); 
unzip(elev_data[0]); 
cp_from_cd(elev_data  [  1  ]); 
unzip(elev_data[l]); 
cp_from_cd(elev_data[2]); 
unzip(elev_data  [2]); 
cp_from_cd(elev_data  [3]); 
unzip(elev_data[3]); 

if(!(fpl= 

fopen(elev_data[0],"r"))){ 

printf("USGS  DEM  file  %s  unavailable  ...  exiting  ...\n",elev_data[0]); 
exit(O); 

} 

if(!(fp2= 

fopen(elev_data[l],"r"))){ 

printf("USGS  DEM  file  %s  unavailable  ...  exiting  ...\n'',elev_data[l]); 
exit(O); 


if(!(fp3= 

fopen(elev_data[2],"r"))){ 

printf("USGS  DEM  file  %s  unavailable  ...  exiting 
exit(O); 


) 


if(!(fp4= 

fopen(elev_data[3],"r"))){ 

printf("USGS  DEM  file  %s  unavailable  ...  exiting 

. . .  \n",elev_data[3]); 

exit(O); 

} 

j_lat=-l; 

j=2‘^jmax-l; 

while(j>=l){ 

if(j>=jmax){ 

for(i=0;i<imax;i++){ 

fread(&hi[i], 

sizeof(signed  short),l,fpl); 
if(hi[i]  ==  -9999) 
hi[i]=0; 

) 

for(i=imax;i<2*imax;i++)  { 
fread(&hi[i], 

sizeoffsigned  short),l/fp2); 
if(hi[i]  ==  -9999) 
hi[i]=0; 


else( 

for(i=0;i<imax;i++){ 

fread(&hi[i], 

sizeoffsigned  short),l,fp3); 
if(hi[i]  ==  -9999) 
hi[i]=0; 

} 

for(i=imax;i<2*imax;i++){ 

fread(&hi[i], 

sizeoffsigned  short),l/fp4); 
if(hi[i]  ==  -9999) 
hi[i]=0; 

} 

} 

get_datal(xminlat,xmaxlat, 
xminlon,xmaxlon,0,2*imax, 
lat_ol  [j],lat_ol  [j-l],&index_r); 
if(index_r  ==  2) 
goto  labell310; 
j--; 


Iabell310: 

fclose(fpl); 

rm_file(elev_data[0]); 

fclose(fp2); 

rm_file(elev_data[l]); 

fclose(fp3); 

rm_file(elev_data[2]); 

fclose(fp4); 

rm_file(elev_datal3]); 

} 

} 

else  if(index_region  ==  2){ 

jmaxl=n2; 

imaxl=nl; 

jmax2=n4; 

imax2=n3; 

/*  1-tile  V 
if(nk  ==  0){ 
jmax=jmax2; 
imax=imax2; 
goto  labelSOOO; 

} 

/*  2-tiles  V 
else  if(nk  ==  1){ 

if(dlat0[0]==dlat0[l]  &&  dlatlfO] 

==dlatl[l]){ 

jmax=jmax2; 

imax=imax2; 

goto  labelbOOO; 

} 

else{ 

jmaxl=n2; 
imaxl=nl; 
jmax2=n4; 
imax2=n3; 
for(j=0;j  <jmaxl  ;j++) 

/*  top  tile  */ 
lat_ol[j]=dlat0[0]+ 
center*(float)j; 
for(j=0;j<jmax2;j++) 

/*  bottom  tile  */ 
Iat_o2[j]=dlat0[l]+ 
center*(float)j; 
for(i=0;i<imaxl;i++) 

/*  top  tile  */ 
lon_ol  [i]=dlon0[0]+ 
center*(float)i; 
for(i=0;i<imax2;i++) 


I*  bottom  tile  * ! 

Ion_o2[i]=dlon0[l]+ 

center*(float)i; 

I*  open  elevation  data  * ! 

cp_from_cd(elev_data[0]); 

unzip(elev_data[0]); 

cp_from_cd(elev_data[l]); 

unzip(elev_data[l]); 

if(!(fpl= 

fopen(elev_data[0],"r"))){ 

printfC'USGS  DEM  file  %s  unavailable  ...  exiting  ...\n",elev_data[0]); 
exit(O); 

} 

if(!(fp2= 

fopen(elev_data[l],"r"))){ 

printffUSGS  DEM  file  %s  unavailable  ...  exiting  ...\n",elev_data[l]); 
exit(O); 

} 

j_lat=-l; 
j=jmaxl-l; 
while(j>=l){ 
f  or(i=0;i<imaxl;i++)  { 
fread(&hi[i], 

sizeoffsigned  short),l/fpl); 

if(hi[i]==-9999) 

hi[i]=0; 

} 

get_datal(xinirJat,xmaxlat, 

xminlon,xmaxlon,0,imaxl, 

lat_ol[j], 

lat_ol  [j-l],&index_r); 
if(index_r  ==  2) 
goto  label3310; 
j-r 

} 

label3310: 

fclose(fpl); 

rm_file(elev_data[0]); 

j=jmax2-l; 

while(j>=l){ 

for(i=0;i<imax2;i++){ 

fread(&hi[i], 

sizeoffsigned  short),l/fp2); 
if(hi[i]  ==  -9999) 
hi[i]=0; 

) 

if(index_r  ==  2) 


goto  label4310; 

j--; 

} 

label4310: 

fclose(fp2); 

rm_file(elev_data[l]); 

} 


/*  3-tiles  V 
else  if(nk  ==  2){ 
jmaxl=n2; 
imaxl=nl; 

jmax2=n4; 

imax2=n3; 

/*  open  elevation  data  * ! 

cp_from_cd(elev_data[0]); 

unzip(elev_data[0]); 

cp_from_cd(elev_data[l]); 

unzip(elev_data[l]); 

cp_from_cd(elev_data[2]); 

unzip(elev_data[2]); 

if(!(fpl=fopen(elev_data[0],"r"))){ 

printf("USGS  DEM  file  %s  unavailable  ...  exiting  ...\n",elev_data[0]); 
exit(O); 

} 

if(!(fp2=fopen(elev_data[l],"r"))){ 

printf("USGS  DEM  file  %s  unavailable  ...  exiting  ...\n",elev_data[l]); 
exit(O); 

} 

if(!(fp3=fopen(elev_data[2],"r"))){ 
printf("USGS  DEM  file  %s  unavailable  ...  exiting 
...\n",elev_data[2]); 
exit(O); 

) 

if  (dlat0[0]==dlat0[l])  { 
for(j=0;j  <jmaxl  ;j++) 

/*  top  tile  * ! 
lat_ol[j]=dlat0[0]+ 

center*(float)j; 

for(j=0;j<jmax2;j++) 

/*  bottom  tile  */ 

Iat_o2[j]=dlat0[2]+ 

center*(float)j; 


for(i=0;i<2*imaxl;i++) 

I*  top  tile  */ 
lon_ol[i]=dlon0[0]+ 
center*(float)i; 
for(i=0;i<imax2;i++) 

I*  bottom  tile  * ! 

Ion_o2[i]=dlon0[2]+ 

center*(float)i; 

jjat=-l; 

j=jmaxl-l; 

while(j>=l){ 

for(i=0;i<imaxl;i++){ 

fread(&hili], 

sizeof (signed  short),l/fpl); 
if(hi[i]  ==  -9999) 
hi[i]=0; 

} 

for(i=imaxl; 

i<2*imaxl;i++){ 

fread(&hi[i], 

sizeof(signed  short),l/fp2); 
if(hi[i]  ==  -9999) 
hi[i]=0; 

} 

get_datal  (xminlat,xmaxlat, 

xminlon,xniaxlon,0,2’^imaxl,lat_ol[j], 

lat_ol[j-l],&index_r); 

if(index_r  ==  2) 

goto  label5310; 

} 


label5310: 

fclose(fpl); 

rm_file(elev_data[0]); 

fclose(fp2); 

rm_file(elev_data[l]); 

j=jmax2-l; 
while(j>=l){ 
f  or(i=0;i<imax2;i+H-)  { 
fread(&hi2[i], 
sizeof(signed  short),l/fp3); 
if(hi2[i]  ==  -9999) 
hi2[i]=0; 

} 

if(index_r  ==  2) 
goto  label6310; 

} 
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Iabel6310: 

fdose(fp3); 

rm_file(elev_data[2]); 

} 

elseif(dlatlll]  ==  dlatl[2]){ 
if(lon_index  ==  1){ 
for(j=0;j<jmaxl;j++) 

I*  top  tile  */ 
lat_ol[j]=dlat0[0]+ 
center*(float)j; 
for(j=0;j<jmax2;j++) 

/*  bottom  tile  */ 

Iat_o2[j]=dlat0[l]+ 

center*(float)j; 

for(i=0;i<imaxl;i++) 

I*  top  tile  * ! 

Ion_ol[i]=dlon0[0] 

+  center  *  (float)i; 
for(i=0;i<2*imax2;i++) 

/*  bottom  tile  *  / 

Ion_o2[i]=dlon0[l] 

+  center  *  (float)i; 

} 

else  if(lon_index  ==  2){ 
for(j=0;j<jmaxl  ;j++) 

I*  top  tile  * ! 
lat_ol[i]=dlat0[0] 

+  center  *  (float)j; 
for(j=0;j<jmax2;j++) 

/*  bottom  tile  *  / 

Iat_o2[j]=dlat0[l] 

+  center  *  (float)j; 
for(i=0;i<imaxl;i++) 

/*  top  tile  * ! 
lon_ol  [i]=dlon0[0] 

+  center  *  (float)i; 
for(i=0;i<2*imax2;i++) 

I*  bottom  tile  * ! 

Ion_o2[i]=dlon00[l]  +  center  *  (float)i; 

) 

jjat=-l; 

j=jmaxl-l; 

while(j>=l){ 

for(i=0;i<imaxl;i++){ 

fread(&hi[i]. 


sizeof(signed  short),l/fpl); 
if(hi[i]  ==  -9999) 
hi[i]=0; 

} 

get_data  1  (xmirLlat,xmaxlat, 
xminlon,xmaxlon,0,2*innax, 
lat_ol[j], 

lat_ol  [j-1  ],&index_r); 
if(index_r  ==  2) 
goto  label7310; 

H 


label7310: 

fclose(fpl); 

rm_file(elev_data[0]); 

j=jmax2-l; 

while(j>=l)| 

for(i=0;i<imax2;i++){ 

fread(&hi[i], 

sizeof(signed  short),l/fp2); 
if(hi[i]  ==  -9999) 
hi[i]  =  0; 

) 

for(i=imax2; 

i<2*imax2;i++){ 

fread(&hi[i], 

sizeof(signed  short),l/fp3); 
if(hi[i]  ==  -9999) 
hi[i]  =  0; 

} 

if(index_r  ==  2) 
goto  label7410; 

} 

label7410: 

fclose(fp2); 

rin_file(elev_data[l]); 

fclose(fp3); 

rm_file(elev_data[2]); 

} 

) 

else  if(nk  ==  3){ 

/*  4-tiles  */ 

cp_from_cd(elev_data[0]); 

unzip(elev_data[0]); 


cp_from_cd(elev_data  [  1  ]); 

unzip(elev_data[l]); 

cp_from_cd(elev_data[2]); 

unzip(elev_data[2]); 

cp_from_cd(elev_data[3]); 

unzip(elev_data[3]); 

if(!(fpl=fopen(elev_data[0],"r"))){ 

printfC'USGS  DEM  file  %s  unavailable  ...  exiting  ...\n",elev_data[0]); 
exit(O); 

} 

if(!(fp2=fopen(elev_data[l],"r"))){ 

printfC'USGS  DEM  file  %s  unavailable  ...  exiting  ...\n",elev_data[l]); 
exit(O); 

) 

if(!(fp3=fopen(elev_data[2],"r"))){ 

printfC'USGS  DEM  file  %s  unavailable  ...  exiting  ...\n",elev_data[2]); 
exit(O); 

) 

if(!(fp4=fopen(elev_data[3],"r"))){ 

printfC'USGS  DEM  file  %s  unavailable  ...  exiting  ...\n",elev_data[3]); 
exit(O); 

) 

for(j=0;j<jmaxl;j++) 
lat_ol[j]=dlatO[0]+ 
center*(float)j;  /*  top  tile  */ 
f  or(j =0;j  <jmax2;j ++) 

Iat_o2[j]=dlat0[2]+ 
center*(float)j;/*bottom  tile  */ 
for(i=0;i<2*imaxl;i++) 
lon_ol[i]=dlon0[0]+ 
center*(float)i;  I*  top  tile  * ! 
for(i=0;i<2*imax2;i++) 

Ion_o2[i]=dlon0[2]+ 
center’^(float)i;/*bottom  tile  * ! 
j_lat=-l; 
j=jmaxl-l; 

■while(j>=l){ 
for(i=0;i<imaxl  ;i++){ 
fread(&hi[i], 

sizeof(signed  short),l/fpl); 
if(hi[i]  ==  -9999) 
hi[i]=0; 

} 

for(i=imaxl;i<2*imaxl;i++){ 

fread(&hi[i], 

sizeoffsigned  short),l/fp2); 
if(hi[i]  ==  -9999) 
hi[i]=0; 


) 

get_datal  (xminlat,xmaxlat, 
xmin]on,xmaxlon,0,2*imaxl, 
lat_ol  [j],lat_ol  [j-l],&index_r); 
if(index_r  ==  2) 
goto  label8310; 

j-; 

} 

label8310: 

fclose(fpl); 

rm_file(elev_data[0]); 

fclose(fp2); 

rm_file(elev_data[l]); 

j=jmax2-l; 

while(j>=l){ 

f  or(i=0;i  <imax2;i++)  { 

fread(&hi2[i],sizeof(signed  short),l/fp3); 

if(hi2[i]  ==  -9999) 

hi2[i]  =  0; 

} 

f  or(i=imax2;i<2*imax2;i++)  { 
fread(&hi2[i],sizeof(signed  short),l/fp4); 
if(hi2[i]  ==  -9999) 
hi2[i]  =  0; 

} 

if(index_r  ==  2) 
goto  label8410; 

j--; 

) 

label8410: 

fclose(fp3); 

rm_file(elev_data[2]); 

fclose(fp4); 

rm_file(elev_data[3]); 

} 

} 

I*  reverse  the  order  of  the  elevation  data 
in  the  y-direction  */ 

while(j  >=  0){ 

jj++; 

for(i=0;i<i_lon;i-t-+){ 

z[i][jj]=hi_2[i][j]; 

lat[jj]=lat_o_2[j]; 
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lon[i]=lon_o_2[i]; 

} 

j--; 

} 

I*  interpolation  of  z  to  model  lat  I  Ion  */ 

zinterp(i_lon,j_lat,lat,lon); 

fpout=fopen("terrain_data","w"); 

fprintf(fpout,"  %11.6f  %9.6f 

%d  %d  %6.1f\n",origlon,lato,ipl,jpl,gdis); 

for(i=0;i<ipl;i++){ 

for(j=0;j<jpl;j++){ 

if(j==jpl-l) 

fprintf(fpout,"  %4d\n",(int) 

(zs[i][j]+0.5)); 

else 

fprintf(fpout,"  %4d",(int) 

(zs[i][j]+0.5)); 

) 

} 

fclose(fpout); 

fpcheck=fopen("terrain_data","r"); 
fscanf(fpcheck,"%f  %f 

%d  %d  %f",&longi,&lati,&numx,&numy,&gridsp); 

min_val=99999; 

max_val=-99999; 

for(i=0;i<ipl;i++){ 

for(j=0;j<jpl;j++){ 

fscanf(^check,"%d",&val); 

if(val<min_val)min_val=val; 

if(val>max_val)max_val=val; 

} 

} 

fclose(fpcheck); 

printf("\ninterpolated  min=%dm 
interpolated  max=%dm\n",min_val,max_val); 

) 

void  zinterpfint  i_l,int  j_l,float  lat[],float  lon[]) 

{ 

int  i,j,m,n; 
float  x,y; 

float  xl,x2,x3,x4,yl,y2,y3,y4,zl,z2,z3,z4; 
float  zd,zdd,yd,ydd; 

int  foxind; 
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float  biggest=-99999.,smallest=99999.; 
float  big=-99999.,small=99999.; 


i=-l; 

j="l; 

for(j=0;j<jpl,-j++){ 

y=lat_b[j]; 

f  or(i=0;i<ipl  ;i++)  { 

x=lon_b[i][j]; 

n=-l; 

do{ 

m=-l; 

do{ 

in++^ 

if(lat_b[j]>= 
lat[n]  &&  lat_b[j]  < 
lat[n+l]  && 
lon_bIi][j]>= 

lon[m]  &&  lon_bIi][}l<lon[m+l]){ 
if  (z[in]  [n]  <smallest) 
smallest=z[m]  [n]; 
if(z[m][n]>biggest) 
biggest=z[m][n]; 

xl=lon[m]; 

yl=lat[n]; 

x2=lon[m+l]; 

y2=lat[n]; 

x3=lon[m]; 

y3=lat[n+l]; 

x4=lon[m+l]; 

y4=lat[n+l]; 

zl=z[in][n]; 

z2=z[m+l]In]; 

z3=z[m][n+l]; 

z4=z[m+l][n+l]; 

yd=yl+((y2-yl) 

/(x2-xl))*(x-xl); 

ydd=y3+((y4-y3) 

/(x4-x3))*(x-x3); 

zd=zl+((z2-zl) 

/(x2-xl))»(x-xl); 

zdd=z3+((z4-z3) 

/(x4-x3))*(x-x3); 


zs[i][j]=zd+ 

(zdd-zd)/  (ydd-yd)*(y-yd); 
found=l; 

} 

else 

found=0; 

}while(m<i_l-l  &&  Ifound); 

)while(n<j_l-l  &&  Ifound); 

} 

} 

} 

void  get_datal (float  xminlat,float  xmaxlat,float  xminlon,float 
xmaxlon,int  io,int  ie,float  lat_l, float  lat_2,int  *index_r) 

{ 

int  i; 

*index_r=l; 

if(lat_l  >=  xmaxlat  &&  lat_2  <  xmaxlat){ 
j_lat+=l; 

lat_o_2[j_lat]=lat_l ; 
i_lon=-l; 

for(i=io;i<ie-l;i++){ 
if(lon_ol[i]  <=  xminlon 
&&  lon_ol[i+l]  >  xminlon)! 
i_lon+=l; 

hi_2[i_lon]lj_lat]=hi[i]; 

lon_o_2[i_lon]=lon_ol[i]; 

} 

else  if(lon_ol[i]  >  xminlon 
&&  lon_ol[i+l]  <  xmaxlon)! 
i_lon+=l; 

hi_2[i_lon][j_lat]=hi[i]; 
lon_o_2[i_lon] =lon_ol  [i]; 

) 

else  if(lon_ol[i]  >=  xmaxlon)! 
i_lon+=l; 

hi_2[i_lon]  [j_lat]=hi[i]; 
lon_o_2[i_lon]=lon_ol  [i]; 
break; 

} 


else  if(lat_l  <  xmaxlat  &&  lat_2  >  xminlat)! 
j_lat+=l; 

lat_o_2  [j_lat] =lat_l ; 
i_lon=-l; 

for(i=io;i<ie-l  ;i++)! 
if(lon_ol[i]  <=  xminlon 


&&  lon_ol[i+l]  >  xminlon){ 
i_lon+=l; 

hi_2[Llon]  [j_lat]  =hi[i]; 
lon_o_2[i_lon]=lon_ol  [i]; 

} 

else  if(lon_ol[i]  >  xminlon 
&&  lon_ol[i+l]  <  xmaxlon){ 
i_lon+=l; 

hi_2[i_lon]  [j_lat] =hi[i]; 
lon_o_2[i_lon]=lon_ol  [i]; 

} 

else  if(lon_ol[i]  >=  xmaxlon){ 
i_lon+=l; 

hi_2[i_lon]  [j_lat]  =hi[i]; 
lon_o_2[i_lon]=lon_ol  [i]; 
break; 

} 

} 

} 

else  if(lat_l  <=  xminlat){ 
j_lat+=l; 

lat_o_2[j_lat]=lat_l; 

i_lon=-l; 

for(i=io;i<ie-l;i++){ 
if(lon_ol[i]  <=  xnainlon  && 
lon_ol[i+l]  >  xmiiilon){ 
i_lon+=l; 

hi_2[i_lon][j_lat]=hi[i]; 
lon_o_2[i_lon]=lon_ol  [i]; 

} 

else  if(lon_ol[i]  > 

xminlon  &&  lon_ol[i+l]  <  xmaxlon){ 
i_lon+=l; 

hi_2[i_lon]  [j_lat]  =hi[i]; 
lon_o_2[i_lon]=lon_ol  [i]; 

} 

else  if(lon_ol[i]  >=  xmaxlon){ 
i_lon+=l; 

hi_2[i_lon][j_lat]=hi[i]; 

lon_o_2[i_lon]=lon_ol[i]; 

*index_r=2; 

break; 


} 

} 

} 
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void  get_data2(float  xminlat,float  xmaxlat,float  xminlon,float 
xmaxlon.int  io4nt  ie,float  lat_l,float  lat_2,int  *index_r) 

*index_r=l; 


if(lat_l  >=  xmaxlat  &&  lat_2  <  xmaxlat)| 
j_lat+=l; 

lat_o_2[j_lat]=lat_l ; 
i_lon=-l; 

f  or(i=io;i  <ie-l  ;i++)  { 
if(lon_o2[i]  <=  xminlon  && 
lon_o2[i+l]  >  xminlon){ 
i_lon+=l; 

hi_2[i_lon][j_lat]=hi2[i]; 

lon_o_2[i_lon]=lon_o2[i]; 

} 

else  if(lon_o2[i]  > 

xminlon  &&  lon_o2[i+l]  <  xmaxlon){ 
i_lon+=l; 

hi_2[i_lon][j_lat]=hi2[i]; 

lon_o_2[i_lon]=lon_o2[i]; 

} 

else  if(lon_o2[i]  >=  xmaxlon){ 
i_lon+=l; 

hi_2[i_lon][j_lat]=hi2[i]; 
lon_o_2  [i_lon] =lon_o2  [i]; 
break; 

} 

) 

} 

else  if(lat_l  <  xmaxlat  &&  lat_2  >  xminlat)i 
j_lat+=l; 

lat_o_2[j_lat]=lat_l ; 
i_lon=-l; 

for(i=io;i<ie-l;i++){ 
if(lon_o2[i]  <=  xminlon  && 
lon_o2[i+l]  >  xminlon)! 
i_lon+=l; 

hi_2[i_lon]  [j_lat]  =hi2[i]; 
lon_o_2[i_lon]=lon_o2[i]; 

) 

else  if(lon_o2[i]  > 

xminlon  &&  lon_o2[i+l]  <  xmaxlon)( 
i_lon+=l; 

hi_2[i_lon][j_lat]=hi2[i]; 

lon_o_2[i_lon]=lon_o2[i]; 

) 


else  if(lon_o2[i]  >=  xmaxlon){ 


i_lon+=l; 

hi_2Ii_lon][j_lat]=hi2[i]; 

lon_o_2[i_lon]=lon_o2[i]; 

break; 

} 


else  if(lat_l  <  xmaxlat  &&  lat_2  >  xirurJat)! 
j_lat+=l; 

lat_o_2[j_lat]=lat_l; 

i_lon=-l; 

for(i=io;i<ie-l;i++){ 
if(lon_o2[i]  <=  xrrdnlon  && 
lon_o2[i+l]  >  xminlon){ 
i_lon+=l; 

hi_2[i_lon][j_lat]=hi2[i]; 

lon_o_2[i_lon]=lon_o2[i]; 

) 

else  if(lon_o2[i]  > 

xminlon  &&  lon_o2[i+l]  <  xmaxlon){ 
i_lon+=l; 

hi_2[i_lon][j_lat]=hi2[i]; 

lon_o_2[i_lon]=lon_o2[i]; 

) 

else  if(lon_o2[i]  >=  xmaxlon){ 
i_lon+=l; 

hi_2[i_lon]  [j_lat] =hi2[i]; 

lon_o_2[i_lon]=lon_o2[i]; 

break; 

} 

} 

} 

else  if(lat_l  <=  xminlat){ 
j_lat+=l; 

lat_o_2[j_lat]=lat_l; 

i_lon=-l; 

f  or(i=io;i<ie-l ;!++){ 
if(lon_o2[i]  <=  xminlon  && 
lon_o2[i+l]  >  xininlon){ 
i_lon+=l; 

hi_2[i_lon]  [j_lat] =hi2[i]; 
lon_o_2[i_lon]=lon_o2[i]; 

} 

else  if(lon_o2[i]  > 

xminlon  &&  lon_o2[i+l]  <  xmaxlon){ 
i_lon+=l; 

hi_2[i_lon][j_lat]=hi2[i]; 

lon_o_2[i_lon]=lon_o2[i]; 
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} 

else  if(lon_o2[i]  >=  xmaxlon){ 
iJon+=l; 

hi_2[i_lon][j_lat]=hi2[i]; 

lon_o_2[i_lon]=lon_o2[i]; 

*index_r=2; 

break; 

) 

) 

) 

} 

void  cp_from_cd(char  str[]) 

{ 

char  system_call[60]; 

strcpy(system_call,"cp "); 

strcat(system_call/Vcdrom/990701_1519/"); 

strcat(system_call,str); 

strcat(system_call,".gz"); 

strcat(system_call/' 

system(system_call) ; 

} 

void  urizip(char  str[]) 


{ 

char  system_calll60]; 
strcpy(system_call,"gzip  -d "); 
strcat(system_call,str) ; 
strcat(system_call,".gz"); 
system(system_call); 

} 

void  rm_file(char  str[]) 


char  system_call[60]; 
strcpy(system_call,"rm  -f "); 
strcat(system_call,str) ; 
system(system_call); 


NASA  MARSHALL  SPACE  FLT  CTR 
ATMOSPHERIC  SCIENCES  DIV 
E501 

ATTN  DR  FICHTL 
HUNTSVILLE  AL  35802 

NASA  SPACE  FLT  CTR 
ATMOSPHERIC  SCIENCES  DIV 
CODE  ED  41  1 
HUNTSVILLE  AL  35812 

US  ARMY  MISSILE  CMND 
AMSMl  RD  AC  AD 
ATTN  DR  PETERSON 
REDSTONE  ARSENAL  AL  35898-5242 

US  ARMY  MISSILE  CMND 
AMSMl  RD  AS  SS 
ATTN  MR  H  F  ANDERSON 
REDSTONE  ARSENAL  AL  35898-5253 

US  ARMY  MISSILE  CMND 
AMSMl  RD  AS  SS 
ATTN  MR  B  WILLIAMS 
REDSTONE  ARSENAL  AL  35898-5253 

US  ARMY  MISSILE  CMND 
AMSMl  RD  DE  SE 
ATTN  MR  GORDON  LILL  JR 
REDSTONE  ARSENAL  AL  35898-5245 

US  ARMY  MISSILE  CMND 
REDSTONE  SCI  INFO  CTR 
AMSMl  RD  CS  R  DOC 
REDSTONE  ARSENAL  AL  35898-5241 

US  ARMY  MISSILE  CMND 
AMSMl 

REDSTONE  ARSENAL  AL  35898-5253 

PACIFIC  MISSILE  TEST  CTR 
GEOPHYSICS  DIV 
ATTN  CODE  3250 
POINT  MUGU  CA  93042-5000 

NAVAL  OCEAN  SYST  CTR 
CODE  54 

ATTN  DR  RICHTER 
SAN  DIEGO  CA  52152-5000 

METEOROLOGIST  IN  CHARGE 
KWAJALEIN  MISSILE  RANGE 
PO  BOX  67 

APO  SAN  FRANCISCO  CA  96555 


DEPT  OF  COMMERCE  CTR 
MOUNTAIN  ADMINISTRATION 
SPPRT  CTR  LIBRARY  R  51 
325  S  BROADWAY 
BOULDER  CO  80303 

DR  HANS  J  LIEBE 
NTIA  ITS  S  3 
325  S  BROADWAY 
BOULDER  CO  80303 

NCAR  LIBRARY  SERIALS 
NATL  CTR  FOR  ATMOS  RSCH 
PO  BOX  3000 
BOULDER  CO  80307-3000 

DEPT  OF  COMMERCE  CTR 
325  S  BROADWAY 
BOULDER  CO  80303 

HEADQUARTERS  DEPT  OF  ARMY 
DAMI  POI 
ATTN  LEE  PAGE 
WASHINGTON  DC  203)0-1067 

MIL  ASST  FOR  ENV  SCI  OFC 
OF  THE  UNDERSEC  OF  DEFNS 
FOR  RSCH  &  ENGR  R&AT  E  LS 
PENTAGON  ROOM  3D129 
WASHINGTON  DC  20301-3080 

DEAN  RMD 
ATTN  DR  GOMEZ 
WASHINGTON  DC  20314 

US  ARMY  INFANTRY 
ATSH  CD  CS  OR 
ATTN  DR  E  DUTOIT 
FT  BENNING  GA  30905-5090 

HQ  AFWA/DNX 

106  PEACEKEEPER  DR  STE  2N3 
OFFUTT  AFB  NE  68113-4039 

PHILLIPS  LABORATORY 
PLLYP 

ATTN  MR  CHISHOLM 
HANSCOM  AFB  MA  01731-5000 

ATMOSPHERIC  SCI  DIV 
GEOPHYISCS  DIRCTRT 
PHILLIPS  LABORATORY 
HANSCOM  AFB  MA  01731-5000 

PHILLIPS  LABORATORY 
PL  LYP  3 

HANSCOM  AFB  MA  01731-5000 
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US  ARMY  MATERIEL  SYST 
ANALYSIS  ACTIVITY 
AMXSY 

ATTN  MR  H  COHEN 
APGMD  21005-5071 

US  ARMY  MATERIEL  SYST 
ANALYSIS  ACTIVITY 
AMXSY  AT 
ATTN  MR  CAMPBELL 
APG  MD  21005-5071 

US  ARMY  MATERIEL  SYST 
ANALYSIS  ACTIVITY 
AMXSY  CR 
ATTN  MR  MARCHET 
APGMD  21005-5071 

ARL  CHEMICAL  BIOLOGY 
NUC  EFFECTS  DIV 
AMSRL  SL  CO 
APGMD  21010-5423 

US  ARMY  MATERIEL  SYST 
ANALYSIS  ACTIVITY 
AMSXY 

APGMD  21005-5071 

ARMY  RESEARCH  LABORATORY 
AMSRL  D 

2800  POWDER  MILL  ROAD 
ADELPHI  MD  20783-1145 

ARMY  RESEARCH  LABORATORY 
AMSRL  OP  Cl  SD  TL 
2800  POWDER  MILL  ROAD 
ADELPHI  MD  20783-1 145 

ARMY  RESEARCH  LABORATORY 
AMSRL  CILL 
ADELPHI  MD  20703-1 197 

ARMY  RESEARCH  LABORATORY 
AMSRL  SS  SH 
ATTN  DR  SZTANKAY 
2800  POWDER  MILL  ROAD 
ADELPHI  MD  20783-1 145 

ARMY  RESEARCH  LABORATORY 
AMSRL  IS 
ATTN  J  GANTT 
2800  POWDER  MILL  ROAD 
ADELPHI  MD  20783-1197 

ARMY  RESEARCH  LABORATORY 
AMSRL  DD 
ATTN  I  ROCCHIO 
2800  POWDER  MILL  ROAD 
ADELPHI  MD  20783 


ARMY  RESEARCH  LABORATORY 
AMSRL 

2800  POWDER  MILL  ROAD 
ADELPHI  MD  20783-1145 

NATIONAL  SECURITY  AGCY  W21 
ATTN  DR  LONGBOTHUM 
9800  SAVAGE  ROAD 
FT  GEORGE  G  MEADE  MD  20755-6000 

US  ARMY  RSRC  OFC 
ATTN  AMXRO  GS  DR  BACH 
PO  BOX  12211 
RTPNC  27009 

DR  JERRY  DAVIS 
NCSU 

PO  BOX  8208 
RALEIGH  NC  27650-8208 

US  ARMY  CECRL 
CECRL  GP 
ATTNDRDETSCH 
HANOVER  NH  03755-1290 

US  ARMY  ARDEC 
SMCAR IMI 1  BLDG  59 
DOVER  NJ  07806-5000 

ARMY  DUGWAY  PROVING  GRD 
STEDP  MT  DA  L  3 
DUGWAY  UT  84022-5000 

ARMY  DUGWAY  PROVING  GRD 
STEDP  MT  M 
ATTN  MR  BOWERS 
DUGWAY  UT  84022-5000 

DEPT  OF  THE  AIR  FORCE 

OL  A  2D  WEATHER  SQUAD  MAC 

HOLLOMAN  AFB  NM  88330-5000 

PL  WE 

KIRTLAND  AFB  NM  871 18-6008 

USAF  ROME  LAB  TECH 
CORRIDOR  W  STE  262  RL  SUL 
26  ELECTR  PKWY  BLD  106 
GRIFFISS  AFBNY  13441-4514 

AFMC  DOW 

WRIGHT  PATTERSON  AFB  OH  45433-5000 

US  ARMY  FIELD  ARTILLERY  SCHOOL 

ATSFTSMTA 

FT  SILL  OK  73503-5600 
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US  ARMY  FOREIGN  SCI  TECH  CTR 
CM 

220  7TH  STREET  NE 
CHARLOTTESVILLE  VA  22448-5000 

NAVAL  SURFACE  WEAPONS  CTR 
CODE  G63 

DAHLGREN  VA  22448-5000 

US  ARMY  OEC  1 

CSTE  EFS 

PARK  CENTER  IV 

4501  FORD  A VE 

ALEXANDRIA  VA  22302-1458 

US  ARMY  CORPS  OF  ENGRS  1 

ENGR  TOPOGRAPHICS  LAB 
ETL  GS  LB 

FT  BELVOIR  VA  22060 

US  ARMY  TOPO  ENGR  CTR  1 

CETEC  ZC  1 

FT  BELVOIR  VA  22060-5546 

SCI  AND  TECHNOLOGY  1 

101  RESEARCH  DRIVE 
HAMPTON  VA  23666-1340 

US  ARMY  NUCLEAR  CML  AGCY  1 

MONA  ZB  BLDG  2073 
SPRINGFIELD  VA  22150-3198 


USATRADOC  1 

ATCD  FA 

FT  MONROE  VA  2365 1-5170 

ATRC  WSS  R  1 

WSMR  NM  88002-5502 

ARMY  RESEARCH  LABORATORY  1 

AMSRL  IS  S 

INFO  SCI  &  TECH  DIR 

WSMR  NM  88002-5501 

ARMY  RESEARCH  LABORATORY  1 

AMSRL  IS  E 

INFO  SCI  &  TECH  DIR 

WSMR  NM  88002-5501 

ARMY  RESEARCH  LABORATORY  1 

^  AMSRL  IS  W 

INFO  SCI  &  TECH  DIR 
WSMR  NM  88002-5501 

^  DTIC  1 

8725  JOHN  J  KINGMAN  RD 
STE  0944 

FT  BELVOIR  V A  22060-62 1 8 
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US  ARMY  MISSILE  CMND 
AMSMl 

REDSTONE  ARSENAL  AL  35898-5243 

US  ARMY  DUGWAY  PROVING  GRD 
STEDP3 

DUGWAY  UT  84022-5000 

USTRADOC 
ATCD  FA 

FT  MONROE  VA  23651-5170 

WSMR  TECH  LIBRARY  BR 
STEWS  IM  IT 
WSMR  NM  88002 

army  research  laboratory 

AMSRL  SL  EW 
ATTN  MR  HEMNI 
INFO  SCI  &  TECH  DIR 
WSMR  NM  88002-5501 

army  RESEARCH  LABORATORY 

AMSRL  SL  EW 

ATTN  MR  KIRBY 

INFO  SCI  &  TECH  DIR 

WSMR  NM  88002-5501 

Record  copy 

TOTAL 


1 


1 


1 


1 


1 


3 

69 
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